#!/usr/bin/env Rscript

library("plyr")
library("foreach")
library("emIRT")
library("pscl")
library("ggplot2")
library("scales")

## #########
## Set Paths
## #########

p2root <- "../"
p2estimates <- paste(p2root, "data/RDATA/estimates/", sep = "")
p2figs <- paste(p2root, "doc/figs/", sep = "")

## ####
## Data
## ####

load(file.path(p2estimates, "ideal/hou.Rdata"))
lOut1 <- lOut

load(file.path(p2estimates, "../bootstrap/hou.Rdata"))
lOut2 <- lOut

## ########
## Analysis
## ########

vSess <- 102:112

idx <- 1


idx <- 11
out1 <- lOut1[[idx]]
out2 <- lOut2[[idx]]

## save(list = c("out1", "out2"), file = "temp.rda")
## load("temp.rda")

mX1 <- out1$output$output$xbar
mXSD1 <- apply(out1$output$output$x,
               c(2, 3),
               sd
               )

shift1 <- mean(mX1)
dilate1 <- sd(mX1)

mX2 <- out2$output$output$means$x
mXSD2 <- sqrt(out2$output$output$varsBS$x)

polarity <- as.numeric(ifelse(cor(mX1, mX2) > 0, 1, -1))

shift2 <- mean(mX2)
dilate2 <- sd(mX2)

shift1
shift2

dilate1
dilate2

newX1 <- (mX1 - shift1) / dilate1
newX2 <- (mX2 - shift2) / dilate2

## #######
## Mapping
## #######

cols1 <- scales::hue_pal(l = 55)(6)
## shapes1 <- shape_pal(solid = TRUE)(6)

shapes1 <- c(1, 2, 5, 16, 7, 8)

pdf(file = paste(p2figs, "intervals_base.pdf", sep = ""),
    width = 10,
    height = 5
    )

par(mfrow = c(1, 2))

plot(mXSD1 / dilate1,
     mXSD2 / dilate2,
     main = "Standard Error Magnitude across Methods",
     xlab = "IDEAL Standard Error",
     ylab = "EM Bootstrap Standard Error",
     cex = .35,
     pch = 16,
     col = alpha("black", 1),
     ylim = c(0, .2),
     xlim = c(0, .4)
     )
abline(a = 0, b = 1, col = alpha("black", .5), lwd = 1.5)

plot(newX1,
     mXSD1 / dilate1,
     xlab = "Ideal Point",
     ylab = "Standard Error",
     main = "Standard Error Magnitude by Ideal Point",
     pch = shapes1[4],
     cex = .5,
     col = alpha(cols1[4], .75),
     ylim = c(0, .4)
     )
points(newX2 * polarity,
       mXSD2 / dilate2,
       pch = shapes1[1],
       cex = .55,
       col = alpha(cols1[1], .75)
       )
text(x = -1,
     y = 0,
     labels = "EM"
     )
text(x = -1,
     y = .11,
     labels = "IDEAL"
     )

dev.off()


################################################################################
